esort <- function(x, sortvar, ...) {
attach(x)
x <- x[with(x,order(sortvar,...)),]
return(x)
detach(x)
}

##data <- read.csv("divergence.statistics.gimli.2.csv")
data <- read.csv("div.stats.rsm1.csv")

pdf("divergence.plot.feb28.pdf",width=6,height=4,pointsize=10)
sw <- c(1e-06,1e-07,1e-08)
pp <- c(0,0.25,0.5,0.75,1)
par(mfrow=c(3,5),mar=c(2,2,0.5,0.5))
for(i in 1:3){
for(ppi in 1:5){
dd <- data[data["a"]==0&data["p"]==pp[ppi]&data["nu"]==1e-04&data["mu"]==sw[i]&data["b"]==0.02&data["d"]==0.016&data["t"]==7300,]
ddmaxa0 <- dd[dd["max_min"]==1,]
ddmina0 <- dd[dd["max_min"]==-1,]
dd <- data[data["a"]==1&data["p"]==pp[ppi]&data["nu"]==1e-04&data["mu"]==sw[i]&data["b"]==0.02&data["d"]==0.016&data["t"]==7300,]
ddmaxa1 <- dd[dd["max_min"]==1,]
ddmina1 <- dd[dd["max_min"]==-1,]
if(ppi==1){
  plot(0,0,xlim=c(0,7300),ylim=c(0,0.3),pty="s",t="n",xaxs="r",yaxs="r",xlab="Time (years)",ylab="",yaxt="n",las=1,xaxt="n")
  axis(2,at=c(0,0.1,0.2,0.3),lab=c("0","0.1","0.2","0.3"))
} else {
  plot(0,0,xlim=c(0,7300),ylim=c(0,0.3),pty="s",t="n",xaxs="r",yaxs="r",xlab="Time (years)",ylab="",yaxt="n",las=1,xaxt="n")
}
axis(1,at=c(0,1825,3650,5475,7300),lab=c(0,5,10,15,20))
##title(paste("p=",pp[ppi]," mu=",sw[i],sep=""))
#points(t(ddmaxa0["day"]),t(ddmaxa0["divergence"]),t="p",pch=1,col="black")
points(t(ddmina0["day"]),rep(0,length(t(ddmina0["day"]))),t="p",pch="|",col="red")
#points(t(ddmaxa1["day"]),t(ddmaxa1["divergence"]),t="p",pch=19,col="black")
points(t(ddmina1["day"]),rep(0,length(t(ddmina1["day"]))),t="p",pch="|",col="black")
for(repi in 1:3){
dd <- data[data["a"]==0&data["p"]==pp[ppi]&data["nu"]==1e-04&data["mu"]==sw[i]&data["b"]==0.02&data["d"]==0.016&data["t"]==7300&data["rep"]==repi,]
dd <- dd[with(dd,order(day)),]
lines(t(dd["day"]),t(dd["divergence"]),col="red")
}
for(repi in 1:3){
dd <- data[data["a"]==1&data["p"]==pp[ppi]&data["nu"]==1e-04&data["mu"]==sw[i]&data["b"]==0.02&data["d"]==0.016&data["t"]==7300&data["rep"]==repi,]
dd <- dd[with(dd,order(day)),]
lines(t(dd["day"]),t(dd["divergence"]),col="black")
}
}
}
dev.off()


## pdf("divergence_plots_b=0.02_d=0.016_nu=1e-04_t=14600.pdf",paper="letter",pointsize=8)
## sw <- c(1e-06, 1e-07,1e-08)
## pp <- c(0,0.25,0.5,0.75,1)
## par(mfrow=c(3,5),mar=c(5,5,3,0.5))
## for(i in 1:3){
## for(ppi in 1:5){
## dd <- data[data["a"]==0&data["p"]==pp[ppi]&data["nu"]==1e-04&data["mu"]==sw[i]&data["b"]==0.02&data["d"]==0.016&data["t"]==14600,]
## ddmaxa0 <- dd[dd["max_min"]==1,]
## ddmina0 <- dd[dd["max_min"]==-1,]
## dd <- data[data["a"]==1&data["p"]==pp[ppi]&data["nu"]==1e-04&data["mu"]==sw[i]&data["b"]==0.02&data["d"]==0.016&data["t"]==14600,]
## ddmaxa1 <- dd[dd["max_min"]==1,]
## ddmina1 <- dd[dd["max_min"]==-1,]
## plot(0,0,xlim=c(0, 14600),ylim=c(0,0.3),pty="s",t="n",xaxs="r",yaxs="r",xlab="Time (years)",ylab="Divergence",las=1,xaxt="n")
## axis(1,at=c(0, 3650, 7300, 10950, 14600),lab=c(0,10,20,30,40))
## title(paste("p=",pp[ppi]," mu=",sw[i],sep=""))
## points(t(ddmaxa0["day"]),t(ddmaxa0["divergence"]),t="p",pch=1,col="black")
## points(t(ddmina0["day"]),t(ddmina0["divergence"]),t="p",pch=1,col="red")
## points(t(ddmaxa1["day"]),t(ddmaxa1["divergence"]),t="p",pch=19,col="black")
## points(t(ddmina1["day"]),t(ddmina1["divergence"]),t="p",pch=19,col="red")
## for(repi in 1:3){
## dd <- data[data["a"]==0&data["p"]==pp[ppi]&data["nu"]==1e-04&data["mu"]==sw[i]&data["b"]==0.02&data["d"]==0.016&data["t"]==14600&data["rep"]==repi,]
## dd <- dd[with(dd,order(day)),]
## lines(t(dd["day"]),t(dd["divergence"]),lty=2)
## }
## for(repi in 1:3){
## dd <- data[data["a"]==1&data["p"]==pp[ppi]&data["nu"]==1e-04&data["mu"]==sw[i]&data["b"]==0.02&data["d"]==0.016&data["t"]==14600&data["rep"]==repi,]
## dd <- dd[with(dd,order(day)),]
## lines(t(dd["day"]),t(dd["divergence"]),lty=1)
## }
## }
## }
## dev.off()




## pdf("divergence_plots_b=0.02_d=0.005_nu=1e-04_t=14600.pdf",paper="letter",pointsize=8)
## sw <- c(1e-06, 1e-07,1e-08)
## pp <- c(0,0.25,0.5,0.75,1)
## par(mfrow=c(3,5),mar=c(5,5,3,0.5))
## for(i in 1:3){
## for(ppi in 1:5){
## dd <- data[data["a"]==0&data["p"]==pp[ppi]&data["nu"]==1e-04&data["mu"]==sw[i]&data["b"]==0.02&data["d"]==0.005&data["t"]==14600,]
## ddmaxa0 <- dd[dd["max_min"]==1,]
## ddmina0 <- dd[dd["max_min"]==-1,]
## dd <- data[data["a"]==1&data["p"]==pp[ppi]&data["nu"]==1e-04&data["mu"]==sw[i]&data["b"]==0.02&data["d"]==0.005&data["t"]==14600,]
## ddmaxa1 <- dd[dd["max_min"]==1,]
## ddmina1 <- dd[dd["max_min"]==-1,]
## plot(0,0,xlim=c(0, 14600),ylim=c(0,0.3),pty="s",t="n",xaxs="r",yaxs="r",xlab="Time (years)",ylab="Divergence",las=1,xaxt="n")
## axis(1,at=c(0, 3650, 7300, 10950, 14600),lab=c(0,10,20,30,40))
## title(paste("p=",pp[ppi]," mu=",sw[i],sep=""))
## points(t(ddmaxa0["day"]),t(ddmaxa0["divergence"]),t="p",pch=1,col="black")
## points(t(ddmina0["day"]),t(ddmina0["divergence"]),t="p",pch=1,col="red")
## points(t(ddmaxa1["day"]),t(ddmaxa1["divergence"]),t="p",pch=19,col="black")
## points(t(ddmina1["day"]),t(ddmina1["divergence"]),t="p",pch=19,col="red")
## for(repi in 1:3){
## dd <- data[data["a"]==0&data["p"]==pp[ppi]&data["nu"]==1e-04&data["mu"]==sw[i]&data["b"]==0.02&data["d"]==0.005&data["t"]==14600&data["rep"]==repi,]
## dd <- dd[with(dd,order(day)),]
## lines(t(dd["day"]),t(dd["divergence"]),lty=2)
## }
## for(repi in 1:3){
## dd <- data[data["a"]==1&data["p"]==pp[ppi]&data["nu"]==1e-04&data["mu"]==sw[i]&data["b"]==0.02&data["d"]==0.005&data["t"]==14600&data["rep"]==repi,]
## dd <- dd[with(dd,order(day)),]
## lines(t(dd["day"]),t(dd["divergence"]),lty=1)
## }
## }
## }
## dev.off()










